Efficient atomic self-interaction correction scheme for non-equilibrium quantum transport 
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Density functional theory calculations of electronic transport based on local exchange and correlation func- 
tional contain self-interaction errors. These originate from the interaction of an electron with the potential 
generated by itself and may be significant in metal-molecule-metal junctions due to the localized nature of 
the molecular orbitals. As a consequence, insulating molecules in weak contact with metallic electrodes erro- 
neously form highly conducting junctions, a failure similar to the inability of local functionals of describing 
Mott-Hubbard insulators. Here we present a fully self-consistent and still computationally undemanding self- 
interaction correction scheme that overcomes these limitations. The method is implemented in the Green's 
function non-equilibrium transport code Smeagol and applied to the prototypical cases of benzene molecules 
sandwiched between gold electrodes. The self-interaction corrected Kohn-Sham highest occupied molecular 
orbital now reproduces closely the negative of the molecular ionization potential and is moved away from the 
gold Fermi energy. This leads to a drastic reduction of the low bias current in much better agreement with 
experiments. 

PACS numbers: 



Electronic transport in molecular devices is currently an 
area of much research interest. Possible applications of such 
devices include high-performance computer components [H, 
chemical sensors |2], disposable electronics, as well as possi- 
ble medical applications |3]. 

Recently it has become possible to construct single 
molecule junctions and measure their transport properties. 
Methods used to construct such devices include mechanically 
controllable break junctions 13,11, Si, STM tips ||6i,|7], litho- 
graphically fabricated nanoelectrodes and colloid solu- 
tions [9]. However, there is much disagreement between the 
results for different experimental methods, with differences 
in the conductance of up to three orders of magnitude being 
reported for the same molecule For this reason 

ab initio quantum transport algorithms are becoming increas- 
ingly important. 

Most computational methods for calculating electronic 
transport involve the combination of scattering theory in the 
form of the non-equilibrium Green's function (NEGF) formal- 
ism lITotl . with an electronic structure method such as density- 
functional theory (DFT] Other schemes include 
time-dependent DFT 11411 or many-body methods 1 15]. How- 
ever, there is also much disagreement between the different 
theoretical methods, as well as between experiments and the- 
ory. In the prototypical case of benzenedi thiol (BDT) attached 
to gold surfaces, the conductance for the most probable con- 
tact geometry calculated with NEGF and DFT is higher than 
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FIG. 1: Schematic diagram of a metal-molecule-metal junction, 
(a) A scattering region is sandwiched between two current/voltage 
probes kept at the chemical potentials fii and fi2- (b) Energy level 
line up and (c) transmission coefficient as a function of energy. 



Hs, which is used to construct the non-equilibrium Green's 
function G{E) 



G{E) = Km [{E + iri) - - Si - S2] 



(1) 



that of any experiments by at least one order of magnitude 
llH, 12, 18, IMl- Still it is desirable to use DFT over many- 



body approaches since it is conceptually simple and computa- 
tionally both undemanding and scalable. 

A schematic illustration of the NEGF method is presented 
in figure [T] A scattering region is sandwiched between two 
cuiTent/voltage probes kept at different chemical potentials 
Ml/2 — Ep ± eV/2 with E-p the Fermi energy, V the bias 
and e the electron charge. This is described by a Hamiltonian 



where S1/2 are the self-energy for the electrodes. G{E) en- 
ters in a self-consistent procedure ifioi 12, 121 for evaluating 
the energy level line up (see Fig[T]3) and hence the two probe- 
current. This is simply the integral between fii and /i2 (the 
bias window) of the transmission coefficient T{E), which is a 
superposition of resonances located in correspondence to the 
molecular energy levels (the highest occupied molecular or- 
bital - HOMO - and the lowest unoccupied molecular orbital 
- LUMO - in Fig. [TJ;). 

The crucial point is that Hs and thus G{E) are constructed 
from the single-particle states calculated by an associated 
electronic structure theory, which needs to meet several re- 
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quirements [20]. First the single particle levels must closely 
resemble the physical removal energies of the system, i.e. the 
molecular levels and Ep should line up correctly. Secondly, 
it should work both at integer and fractional occupation, en- 
suring correct charging of the molecular levels. Finally the 
matrix elements between the leads and the molecule should 
be calculated accurately. 

DFT in its Kohn-Sham (KS) formulation Ii21il is the most 
widely used electronic structure theory for transport, although 
in principle it does not meet any of the previous require- 
ments. In fact DFT observables are integral quantities of the 
KS eigenvalues (the total energy) or of the KS wave-functions 
(the charge density), with the individual KS orbitals having 
little meaning. There is however an important exception since 
the energy of the KS HOMO (eHOivr n) i s the negative of the 



ionization potential / of the system 11221 12311 . This suggests 
that, at least for moderate bias and HOMO transport, KS the- 
ory can be used for transport calculations 12411 . Unfortunately 
for standard local functionals, such as the local density ap- 
proximation (LDA), enoMO is nowhere near — /. 

Most of the failures of LDA can be traced back to the self- 
interaction (SI), i.e. the spurious interaction of an electron 
with the Hartree and exchange-correlation (XC) potentials 
generated by itself ll25ll . In the case of an exact exchange the- 
ory (for instance Hartree-Fock) the self-Hartree cancels with 
the self-XC potential, however for LDA the cancellation is in- 
complete. This results in the eigenstates of a molecule being 
too high in energy, which translates in erroneously positioned 
peaks in the transmission coefficient. Self-interaction correc- 
tions (SIC) are possible | |25|] . however their implementation in 
typical solid state codes is cumbersome, since the theory be- 
come orbital dependent and the energy minimization cannot 
follow the standard KS scheme. 

These problems become even more serious when SIC is 
combined with the NEGF method since the KS orbitals are 
never individually available. Recently we have implemented 
lE6ll an alternative and approximated method for dealing with 
SIC. This is based on an atomic approximation (ASIC) i EtIi . 
where the SIC orbitals that minimize the energy are assumed 
to be atomic-like. The correction thus becomes atomic and 
no information are needed other than the charge density and 
the ASIC projectors [26]. Importantly ASIC produces single- 
particle energy levels rather close to the experimental molec- 
ular removal energies, and in particular ehomo ~ ™d 
the chomo for negatively charged molecules is close to the 
molecular affinity. As an example ASIC places chomo for 
1,2-BDT at 8.47 eV to compare with the LDA value of 
4.89 eV and the experimental ionization potential ^^8.5 eV 



We have numerically implemented the ASIC method 02611 
in the localized atomic orbital code SIESTA ll29ll . which is 
the DFT platform for our transport code Smeagol lll2ll . and 
carried out calculations for the prototypical Au/Benzene/Au 
molecular devices. We use a double zeta polarized basis set 
fell for carbon and sulfur s and p orbitals, double zeta for 
the Is orbitals of hydrogen and 6s-only double zeta for gold. 



The mesh cut-off is 200 Ry and we consider 500 real and 80 
complex energy points for integrating the Green's function. 

First we consider a 1,4-BDT molecule sandwiched between 
two fee (111) gold surfaces and attached to the gold hollow 
site (figure |2]l. We optimize the distance between the sulfur 
atom and the gold surface to a value of 1.9A, which agrees 
with previous calculations 1 17, 30ll . 




FIG. 2: BDT molecule attached to the hollow site of the Au (ill) 
surface. The sulfur- surface distance is 1.9A. Color code: Au=yellow, 
C=black;, S=dark yellow, H=blue. 

The orbital resolved density of states (DOS), transmission 
coefficient and I-V curves are presented in figure [3] for both 
LDA and ASIC. From the DOS (panels (a) and (b)) it is clear 
that the effect of ASIC is that of shifting the occupied or- 
bitals downwards relative to E-p of the gold. The HOMO- 
LUMO gap is considerably larger than that of the LDA, and 
most importantly in the case of ASIC there is little DOS orig- 
inating from the molecule at Ep. This has profound effects 
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FIG. 3: Transport properties of a BDT molecule attached to the gold 
(111) hollow site. The left plots correspond to LDA and the right 
ones to ASIC. The upper panels are the DOS of the S and C n or- 
bitals ((a) and (b)), the middle are the transmission coefficients as a 
function of energy for various bias ((c) and (d)) and the lower are 
the I-V curves. Figure (f) is a zoom of (e) and compares our results 
with experiments from reference 01 ■ The vertical lines in (c) and (d) 
mark the bias window. 

over the electron transmission. The LDA peaks of T{E) aris- 
ing from occupied orbitals are shifted downwards in energy 



and away from Ep. At variance from LDA (Fig. where 
T{Ep) is dominated by a resonance at ehomo. the ASIC 
transmission (Fig. [3]i) is through the BDT gap and therefore 
it is tunneling-Hke. This resuhs in a drastic reduction of the 
low-bias current when going from LDA to ASIC (Fig. [3^). 
The ASIC-calculated conductance at zero bias is now about 
O.O6G0 (Go = 2e'^/h), compared to 0.23Go of LDA. A con- 
ductance of O.O6G0 is much closer to the value of O.OllGo 
obtained by Xiao et. al. 01 and is actually lower than values 
0.09-0. 14Go obtained by Tsutsuiet. al. |5]. 

Other anchoring configurations to the (111) surface were 
investigated and their zero-bias conductance are shown in Ta- 
ble H] Note that ASIC returns values in the region of about 
O.O6G0 for several different anchoring geometries. Such sta- 
bility with variation of the anchoring structure is important 
as the peaks in the experimental conductance histograms are 
relatively sharp IlTJ. This indicates that the different metal- 
molecule junctions have similar conductances, despite the fact 
that the anchoring of the molecule to the surface/tip may vary. 
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FIG. 4: BDMT molecule attached to the hollow sites of the gold 
(111) surface. The contact geometry is identical to that for BDT. 
Color code: Au=yellow, C=black, S=dark yellow, H=blue. 
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TABLE I: Zero-bias conductance for different configurations of BDT 
on gold (111). Experimental values include 0.01 IGo 1 7] and 0.09Go 
f6]. The anchoring configurations investigated are: hollow site (Ho), 
bridge site (Br), Au adatom (Ad). Ho (30°) describes a hollow site 
with BDT at a 30° angle with respect to the transport direction, and 
the two last rows correspond to asymmetric anchoring to the two 
electrodes, d is distance between the S atom of the thiol group and 
the plane of the gold surface (or to the adatom for Ad). 



FIG. 5: Transport properties of a BDMT molecule attached to the 
gold (111) hollow site. The left panels correspond to LDA and the 
right ones to ASIC. The upper panels are the DOS of the S and C tt 
orbitals ((a) and (b)), the middle are the transmission coefficients as 
a function of energy for various bias ((c) and (d)) and the lower are 
the I-V curves. Figure (f) is a zoom of (e) and compares our results 
with experiments from reference The vertical lines in (c) and (d) 
mark the bias window. 



As a second case we investigate benzenedimethanethiol 
(BDMT) molecules on the same gold surface (Fig. [Hi. Also 
for this molecule we consider a number of anchoring struc- 
tures, although here we report only for the hollow site since 
different structures do not present qualitatively different re- 
sults. 

In contrast to BDT, the DOS of BDMT already presents 
a large HOMO-LUMO gap in LDA (see Fig. |5}i), which is 
further increased by the ASIC This time ASIC has the 
only effect of changing the alignment of enoMO with respect 
to i?F, which cuts through the HOMO-LUMO gap for both 
LDA and ASIC. Therefore both LDA and ASIC offer a pic- 
ture of tunneling-like transport through the molecular gap, and 
little difference can be found (see Fig. [Sj;, d, e and f). Hence, 
the I-V curves are quite similar, with currents approximately 
one order of magnitude smaller than those of BDT. 

With these results in hand we can conclude that in general 



ASIC drastically improves the agreement between theory and 
experiments. However some disagreement still remains. In 
particular it appears that even in the case of ASIC the current 
at low bias is larger than that typically measured. Although 
an exhaustive comparison is complicated by the fact that the 
experimental spread of the data is large, here we speculate on 
the possible source of such disagreement. First one may argue 
that the contact geometry is not correct. Indeed recent X-ray 
standing wave experiments 1I31I1 demonstrate that S atoms in 
thiol groups on gold join more favorably to adatoms. This 
means that the Au-S-molecule moiety may be the one rele- 
vant for the transport experiments. However calculations with 
two Au adatoms as anchoring sites lead to strong pinning of 
EHOMO to the gold Ey il7lll9tl . ASIC does not change this 
feature and the conductance remains large. An asymmetric 
anchoring configuration with one hollow site and one adatom 
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(see last column of table gives us a ASIC conductance of 
0.03 Go, still higher than experiments. Note that the presence 
of the mobile Au-S-molecule moieties seems also difficult to 
conciliate with the relative robustness of the peaks in the con- 
ductance histograms of breaking-junctions |7]. 

Hence we perhaps have to accept the fact that the disagree- 
ment between theory and experiments persists. Notably the 
problem is now that of calculating accurately tunneling ma- 
trix elements between the S and the surface. The problem is 
thus critically dependent on the quality of the wave-function 
and in turn of the actual scattering potential, for which ASIC 
does not offer substantial improvement over the LDA. In par- 
ticular ASIC still overestimates the polarizability of molecules 
II32II . with a quantitatively incorrect prediction of the response 
exchange and correlation field. In addition it is important to 
remark that we have applied the ASIC only to the molecular 
degrees of freedom, without correcting the Au atoms. It is 
thus likely that the Au 6 s orbitals at the surface are too ex- 
tended, leading to a larger current. 

In conclusion, we have demonstrated that a simple SIC 
scheme is able of lowering the energy levels of the occupied 
molecular orbitals, which now resemble closely the actual 
vertical removal energies. This has profound consequences 
over the transport properties of metal/molecule/metal junc- 
tions since spurious resonances at the Fermi level can be re- 
moved, leading to tunneling transport for molecules for which 
LDA erroneously predicts metallic conductance. The agree- 
ment with experiments is thus greatly improved. 
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